%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function evaluates a 1-D spline which appximating the function on
% the node points outside the approximation interval in a linear way
% satisfying the zero-order and first-order continuity at the boundary
% points.
% Author:   Xing Guo (xingguo@umich.edu)
% Date of this version: 2017/9/18
% Function: 
% [B,x]=ExtendedSplineEvaluation(breaks,PP,x,order)
% Input variables: 
%   breaks:     vector of break points
%   PP:         pp form of spline accepted by Matlab function ppval
%   x:          Evaluation points
%   order:      A scaler or a vector for the evaluation order
% Output variables:
%   SpValue:    Evaluation Results (cell if order is a vector)
% Special Note: 
% Uses:     ppval
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [SpValue]=ExtendedSplineEvaluation(breaks,PP,x,order)

SpValue         =   cell(length(order),1);
x_left          =   breaks(1);
x_right         =   breaks(end);
I_left          =   x<x_left;
I_right         =   x>x_right;

for i=1:length(order)
    SpValue{i}      =   zeros(length(x),1);
    if any(I_left)
        y_left      =   ppval(PP.d0,x_left);
        slope_left  =   ppval(PP.d1,x_left);
        switch order(i)
            case 0
                SpValue{i}(I_left)      =   y_left+slope_left*(x(I_left)-x_left);                
            case 1
                SpValue{i}(I_left)      =   slope_left;
        end
    end

    if any(I_right)
        y_right     =   ppval(PP.d0,x_right);
        slope_right =   ppval(PP.d1,x_right);
        switch order(i)
            case 0
                SpValue{i}(I_right)     =   y_right+slope_right*(x(I_right)-x_right);                
            case 1
                SpValue{i}(I_right)     =   slope_right;
        end        
    end
    if any(~I_left & ~I_right)
        d_temp      =   ['d',num2str(order(i))];
        SpValue{i}(~I_left & ~I_right)  =   ppval(PP.(d_temp),x(~I_left & ~I_right));
    end
end

if length(order)==1
    SpValue         =   SpValue{1};
end